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The multi-reference Coupled Cluster method first proposed by Meller et al {J. Chem. Phys. 1996^ has 
been implemented and tested. Guess values of the amplitudes of the single and double excitations (the 
T operator) on the top of the references are extracted from the knowledge of the coefficients of the Multi 
Reference Singles and Doubles Configuration Interaction (MRSDCI) matrix. The multiple parentage problem 
is solved by scaling these amplitudes on the interaction between the references and the Singles and Doubles. 
Then one proceeds to a dressing of the MRSDCI matrix under the effect of the Triples and Quadruples, the 
coefficients of which are estimated from the action of T^. This dressing follows the logics of the intermediate 
effective Hamiltonian formalism. The dressed MRSDCI matrix is diagonalized and the process is iterated to 
convergence. The method is tested on a series of benchmark systems from Complete Active Spaces (CAS) 
involving 2 or 4 active electrons up to bond breakings. The comparison with Full Configuration Interaction 
(FCI) results shows that the errors are of the order of a few milli-hartree, five times smaller than those of the 
CASSDCI. The method is totally uncontracted, parallelizable, and extremely flexible since it may be applied 
to selected MR and/or selected SDCI. Some potential generalizations are briefly discussed. 


I. INTRODUCTION 

In the domain of molecular physics and quantum chem¬ 
istry the many-body problem is perfectly clear as long as 
it is formulated from a single reference. The perturba¬ 
tive expansion of the wave operator and its diagrammatic 
transcription offer a guide to understand the relations 
between the multiplicative structure of the wave func¬ 
tion and the additive structure of the energy. The linked 
cluster theorenP clarifies the questions of the size consis¬ 
tency and of the strict separability into fragments. The 
defects of truncated Configuration Interaction (Cl) are 
well understood and algorith ms h ave been proposed to re¬ 
spect approximately (CEPA-CP^, CEPA-tPEI) ^ or strictly 
((SC)^CI)P the cancellation of unlinked diagrams. The 
Coupled Cluster (CC) methocP®i^ is definitely the most 
elegant formalism and can be considered as the stan¬ 
dard treatment in its CCSD version, or in the CCSD(T) 
versiorfi^ which incorporates the fourth order effect of the 
triply excited determinants. But all these approaches fail 
when one cannot expect that a single determinant will 
represent a reliable starting point to conveniently gener¬ 
ate the wave function. 

This is precisely the situation in many domains. The 
excited states present an intrinsic multi-determinantal 
character, and frequently a multi-configurational char¬ 
acter. So are the magnetic systems in their low en¬ 
ergy states, and the treatment of chemical reactions, in 
which chemical bonds are broken, also requires to con¬ 
sider geometries for which a single determinant picture 
is not relevant. A generalized linked cluster theorem 
has been established by Brandov^^, which gives a con- 
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ceptual guide, but the conditions that must be fulfilled 
for its demonstration (Complete Active Space (CAS) as 
reference space, mono-electronic zero-order Hamiltonian) 
would lead to strongly divergent behaviors of the corre¬ 
sponding perturbative expansion in any realistic molec¬ 
ular problem. Practical computational tools have been 
proposed, most of them being state-specific. One may 
quote second order perturbation expansions based on de¬ 
terminants from selected references (CIPSI)P^*i^, inter¬ 
mediate Hamiltonian dressin^i^, in the so-called shifted 
Bk techniqu^i^. These methods are not strictly size- 
consistent, the conditions to satisfy the strict separability 
of determinant-based expansions require to define sophis¬ 
ticated zero-order HamiltonianJi^. Contracted pertur¬ 
bative expansions, which perturb the multideterminant 
zero-order wave function under the effect of linear com¬ 
binations of outer-space determinants have al so be en pro¬ 
posed. One may quote the CASPT2 methocP^lSO]^ which 
uses a monoelectronic zero-order Hamiltonian, faces in¬ 
truder state problems and is not size consistent, the 
NEVPT2 methocPiHlll which uses a bi-electronic zero- 
order Hamiltonian (the Dyall’s on^^ and is size con¬ 
sistent and intruder-state free, and the method from 
Werneil^, as well as the perturbation derived by Mukher- 
jee et aPSHSai fj-om their MRCC formalism. 

If the CASSCF wave function is considered as the 
counterpart of the single determinant reference the 
CASSDCI is the counterpart of the SDCI, with the same 
size-inconsistence defect, and the research of MRCEPA 
and MRCC has been the subject of intense method¬ 
ological researches for about 20 years, without evident 
success. The cancellation of all unlinked terms in the 
MR expansion (i.e. a MRCEPA or MR(SC)^CI) for¬ 
malism is not an easy taslP^M!^ jf one lets aside the 
MRCC methods that attribute a specific role to a single 
referenc^^, a few state-specific strictly multi-reference 
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CC methods have been proposed, one by one of the 
authors and collaboratorJi^, another one by Mukher- 
jee an d cow orkerj^, a third one in a Brillouin Wigner 
contexlPSlMI We return here on the first proposal which 
had only be tested on a single problem. We shall present 
briefly the method in section 2, then the principle of its 
implementation (section 3), followed by numerical illus¬ 
trations of its accuracy (section 4). The last section will 
discuss the advantages of this formalism, its flexibility 
and possible extensions. 

II. METHOD 
A. Principle 

Let us call |/) the reference determinants, the number 
of which will be called N. The reference space may be 
a CAS, but this is only compulsory if one wants to sat¬ 
isfy the strict-separability property. If not the method 
is applicable to incomplete model spaces as well. The 
projector on the model space is 

Po = J2\^){i\. ( 1 ) 

I 

Let us consider a zero-order wave function restricted to 
the model space, 

( 2 ) 

/ 

This function may be either the eigenfunction of PqHPq^ 

PoHPom=E^\^^), (3) 

or the projection of the eigenvector of the CASSDCI on 
the model space 

\K) = Po\'^cassuci)- (4) 

The CASSDCI wave function is written as 


B. The multi-parentage problem and the extraction of 
guess values of the excitation amplitudes from the 
CASSDCI eigenvector 

One may easily recognize that there exist some degrees 
of freedom in the determination of the wave operators. In 
the single reference CCSD expansion one searches for the 
amplitudes of the excitations sending from the reference 
$0 to the singly and doubly excited determinants. One 
evaluates the amplitudes of the Triples and Quadruples 
as given by the action of on $o, and the eigenequation 
is projected on each of the Singles and Doubles. If the 
number of Singles and Doubles is n, one may write a set 
of n coupled quadratic equations on the amplitudes. But 
it may be more convenient to guess a first evaluation of 
these amplitudes from the coefficients of the Singles and 
Doubles in the SDCI matrix, which may be done in a 
unique manner. From these amplitudes on may obtain 
a guess of the coefficients of the Triples and Quadruples 
and it is convenient (ensuring for instance a better con¬ 
vergence than solving coupled biquadratic equation) to 
write the process as an iterative dressing of the SDCI 
matrix, in the^irit of Intermediate Effective Hamilto¬ 
nian formalisnPSl 

In the multireference context one faces a genealogical 
problem, sometimes called the multiple-parentage prob¬ 
lem. Actually for a state-specific formalism, one has only 
one coefficient for each of the singly and doubly excited 
determinants |i). In principle one may decide that this 
determinant is obtained from each of the references and 
one would write then 

cT = T.dncT (9) 

/ 

but one must find a criterion to define the N dn am¬ 
plitudes from the knowledge of a single coefficient. Re¬ 
turning to a perturbative estimate of the coefficients of 
the Singles and Doubles starting from 'k™, the first-order 
expression of these coefficients 


|4^CASSDCl)=E""l^) + E^rN) (5) 

I i 

where \i) are the Singles and Doubles (the determinants 
of the CASSDCI space which do not belong to the refer¬ 
ence space). We want to follow a Jeziorski-MonkhorslP^ 
expression of the wave operator Cl which is supposed to 
send from the zero-order wave function to the exact one 

Dl^-™) = l^-™) (6) 

as a sum of reference-dependent operators 

ciPo = j2Cij\i){i\. ( 7 ) 

1 

Each of the ^/’s will take an exponential form 

D/ = exp(T/), (8) 


(1) mPi\t) _Y^m m\i) 


( 10 ) 


suggests that the amplitudes of the excitation operators 
from the references to the Singles and Doubles might sat¬ 
isfy 

This scaling had been proposed in ref.®. This condition 
may be expressed as 

dT. = XT{im) (12) 


where the quantity A™ is the inverse of an energy. Re¬ 
injecting this expression in Eq. § leads to 


and each operator Tj will be truncated to the single and 
double excitations, as one does in the CCSD formalism. 


cr = ArEcr(4|i?i*), (is) 

I 
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which defines A™ as 

m 

ym _ _ 

* “ {nm) 


where {{k,l) S (J —> a)} denotes the couples {k,l) for 
which TikTii = ±T/q,. The sign is governed by the per- 
(14) mntation logics. Then one might write the coefficient c™ 
as 


These are the key equations which define guess values 
of the amplitudes of the excitations leading from the 
references to the Singles and Doubles. Notice that we 
only consider amplitndes for the excitations which corre¬ 
spond to physical interactions, and since H is at most 
bi-electronic, one only introduces single- and double¬ 
excitation operators. Finally, we can re-express the 
CASSDCI wave function as 


c = (21) 

I 

D. Dressing of the CASSDCI matrix 

If one considers the eigenequation relative to (i| 


I'f'cASSDCl) — E E 1-^) (15) 

where 

fn\I) = \i). (16) 

C. Coefficients of the Triples and Quadruples 

One may then generate the Triples and Quadruples |a). 
Among them only those which interact with the Singles 
and Doubles (i.e. which are generated by the action of 
H on the Singles and Doubles and which do not belong 
to the CASSDCI space) have to be considered. One may 
find the references with which they present either 3 or 4 
differences in the occupation numbers of the molecular 
orbitals (MOs). These reference determinants may be 
called the grand-parents of |a). The comparison between 
I a) and each of its grand-parents |/) defines the excita¬ 
tion operator from |/) to |a) as a triple or quadruple 
excitation 

fi^\I) = |a) (17) 


((*iij|*) - A-) cT+Y,m\i)cr + Y.m\j)cf 

I 3i^i 

+ E(*l^l«)C=0 (22) 

Ol 

one may decompose the last term 

OL a I 

= (23) 

Introducing the quantities 


{i\A^\I) =Y,dZm\a), 


(24) 


one may write the eigenequation (22) as 

{z\H\z) - E^) cr + E ((*1^1^) + 

/ 

+ E(*l^bK = o 


(25) 


which may be expressed in second quantization as the 
product of 4 (or 3) creation operators and 4 (or 3) anni¬ 
hilation operators 

fia = a\ala^^a)^aeafagah (18) 

The creations run on active and virtual MOs, the an¬ 
nihilations run on active and inactive occupied MOs but 
the number of inactive indices among the creation and/or 
among the particles must be equal to 3 or 4, otherwise 
the determinant would belong to the CASSDCI space. 
Knowing the operator, it may be factorized as the prod¬ 
uct of two complementary double (or single) excitation 
operators in all possible manners (each double excitation 
keeping untouched the Mg value) 

Tia = ^TjkTii = ±fiyTiu = • • • (19) 

Then we may write the contribution to the coefficient of 
I a) issued from the reference |/) as 

dZ = E 


which suggests to treat the effect of the Triples and 
Quadruples as a column dressing of the CASSDCI ma¬ 
trix. A similar idea has been exploited in the single¬ 
reference CCSD context, which may be presented and 
managed as an iterative dressing of the column be¬ 
tween the reference and the Singles and Double^^. The 
Coupled Cluster dressed CASSDCI Hamiltonian may be 

written as A^assdci (^H Dcassdci, which is non- 

Hermitian. Defining the projector on the Singles and 
Doubles as 

-PsD = A::assdci - Po, (26) 


A^assdci^^-Pcassdci = PsbA"^Po (27) 

one may define an equivalent Hermitian dressing A™' in 
the case where one considers the Hermitization of the 
dressed CASSDCI matrix to be desirable, 

(i|A“'|/) = = (f|A'"|/) 


( 20 ) 


(28) 
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provided that one introduces a diagonal dressing of the 
CASSDCI matrix 


(/|A^ 


'\I) = - — 

L- r 




(/|A™'|z)cr 


The 


ces 


diagonalization 


A:assdci ( 


iJ + A*^ 


of the 
-PcASSDCI 
will 


A^assdci + A™'^ -Roassdci will give the 
desired eigenenergy and eigenvector 


(29) 

matri- 

and 

same 


'PcASSDCI ( 
-PcASSDCI ( 


i? + A^ 


i? + A^ 


While generating the parents of |a), its interactions with 
them, {k\H\a), are stored. At this step, the reference 
grand-parents |/) are identified as having 3 or 4 differ¬ 
ences with I a). Then, the excitation operator leading 
from |J) to I a) is expressed in all possible manners as 
products of two complementary single or double excita¬ 
tions. For each couple of complementary excitations, the 
product of the amplitudes is accumulated to compute 
dfa according to equation (20). Finally, the product of 
{k\H\a) with dyj, is accumulated in (fc|A™|/) for each 
parent \k) of |a) according to equation (24). Once the 
loop on the \i)'s is done, all the |a)’s have been gener- 

PcassdciI^Cc) = -E’ccA^assdciI^'cc)^^®^’ column dressing is completed. Then, in 

order to fit with a symmetric diagonalization technique, 
^ one symmetrizes the dressing as mentioned in the pre¬ 

ceding section (equations (28) and (29)). The dressed 
6 \ _ A i.T.m \CASSDCI matrix is diagonalized and the process is re- 

CASSDCil cc) CC CASSDCI I CC/pgg^^g^j convergence of the calculated dressed en¬ 

ergy. From the computational point of view, this process 
requires the storing of the dressing columns which scales 
as N X (riocci^virt)^ where N is the number of determi¬ 
nants in the reference, Uocc and Uvirt are respectively the 
number of occupied an virtual MOs. This amount of 
memory is reasonable, and does not represent a bottle¬ 
neck for the present applications. Regarding the CPU 
time, the costly part concerns the handling of the |a)’s 
which scales as (nocci^virt)"^- Nevertheless, the process is 
perfectly parallelizable as all the work done with the |q!)’s 
generated from \i) does not depend on the other \i)’s. 


(31) 

Of course the process has to be iterated, the resulting 
eigenvector defines new coefficients on both the references 
and the Singles and Doubles, which lead to new ampli¬ 
tudes, new evaluations of the coefficients of the Triples 
and Quadruples, new dressings. Since the eigenvectors of 
the dressed matrices are identical, the two formulations, 
Hermitian or non-Hermitian, converge to the same solu¬ 
tion. The converged solutions are the MRCCSD energy 
and the MRCCSD amplitudes, which define the expo¬ 
nential wave operator. 


III. IMPLEMENTATION 


B. Practical issues 


The proposed algorithm was implemented in the Quan¬ 
tum Packag^^, an open-source series of programs devel¬ 
oped in our laboratory. The bottleneck of this algorithm 
is the determinant comparisons needed to determine the 
excitation operators and phases during the reconstruc¬ 
tion of the genealogy of the |a)’s. This was made pos¬ 
sible thanks to a very efficient implementation of Slater- 
Condon’s rule^. 


A. General structure 


The definition of A™ can lead to numerical instabilities 
when is small. Nevertheless, in such cases the 

contribution of |i) to the post-CAS correlation energy is 
also small, suggesting that one might use a perturbative 
estimate of A™. In practice, we use the perturbative A^ 
according to two different criteria. The first one concerns 
the ratio of the variational coefficient c™ (obtained at a 
given iteration) over its perturbative estimate (see (ll0|)). 

If ^ [0,0.5] then the amplitudes involving |i) are 
determined using the perturbative defined as 


At each iteration step, one first assigns the values of 
the A™ parameters obtained from the eigenvector of the 


(dressed) CASSDCI matrix according to equation (14). 
From these parameters A(” the amplitudes of the single 
and double excitations are uniquely defined. Then one 
loops on the Singles and Doubles |i). On each of them one 
reapplies the excitation operators to generate the |a)’s. 
Those which belong to the CASSDCI space are elimi¬ 
nated. The parents of |a) (that are all the Singles and 
Doubles |A:)’s such that {a\H\k) ^ 0) are generated. If 
one of the |fc)’s has already been considered in the loop 
on the |f)’s (fc < i) , this jo) has been already generated 
and taken into account and must not be double counted. 


A' 


m(pert) 




(32) 


In such situations the coefficient c]" is not determined 
by its interaction with the reference determinants, but 
comes from higher-order effects. The second criterion 
concerns the absolute value of each of the dfi defined ac¬ 


cording to (12). If any of these terms calculated with the 
A™ obtained from the variational calculation (see 

is larger than 0.5, the perturbative is used to 

determine the amplitudes defined as 




2(pert) 


= HuX] 


'i(pert) 


(33) 
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and the working amplitudes dJl are set to . This 

condition avoids numerical instabilities occurring when 
both c™ and are small, and allows us the con¬ 

trol of the maximum value of the amplitudes. As soon 
as along the iterations one of the |i)’s fulfills one of these 
criteria, it will be treated perturbatively in the following 
iterations. This precaution avoids signihcant oscillations 
due to back and forth movements from perturbative to 
variational treatment of the A™. The numerically ob¬ 
served residual oscillations are of the order of magnitude 
of 10“® Ah, which may certainly be attributed to the 
non linear character of the numerical algorithm. Never¬ 
theless, the order of magnitude of the residual oscillations 
is much smaller than the chemical and even spectroscopic 
accuracy. 

IV. NUMERICAL TEST STUDIES 

We decided to test the accuracy and robustness of the 
method on a series of benchmarks, some of which have 
been used in the evaluations of other MRCC proposals 
and of alternative MR approaches. They essentially con¬ 
cern model problems, especially bond breaking problems 
or the treatment of degenerate situations. They require 
to use a CAS with either two electrons in two MOs or 
four electrons in four MOs. In all cases the method con¬ 
verged in a few iterations. A systematic comparison is 
made with FCI estimates, either taken from the literature 
or obtaine d fro m a CIPSI type variation-|-perturbation 
calculatiorfSHsl where the perturbative residue is about 
-6 mAh. Of course the CASSDCI is already a rather 
sophisticated treatment, which takes into account, al¬ 
though in a size-inconsistent manner, the leading correla¬ 
tion effects, both the non-dynamical part in the CAS and 
the dynamical part in the SDCI step. One may expect 
that the improvement brought by the MRCC treatment 
will be significant when the number of important inactive 
double excitations is large. 

In order to have a global view of the performance of the 
here-proposed algorithm, we report for each calculation 
(except the symmetric dissociation of the water molecule) 
potential energy curves, the error to FCI estimate of 
our MRCCSD algorithm together with the CASSDCI. 
Tables showing the error with respect to the FCI esti¬ 
mate of the MRCCSD and CASSDCI are also presented, 
complemented by the total energies of the FCI estimate. 
The non-parallelism error (NPE) is here calculated as 
the difference between the minimum and maximum er¬ 
ror to the FCI estimate. The spectroscopic constants are 
obtained from an accurate fit of the obtained potential 
energy curves with a generalized Morse potential repre¬ 
sentation. The spectroscopic constants reported here are 
the equilibrium distance i?oq in the frequency fceq in 

E\y /and the atomization energy Df. in kcal/mol. 

All the calculations were performed with the Quantum 
Packag^^, an open-source series of programs developed 
in our group. 


A. Single-bond breakings 

The treatment of the breaking of a single bond in prin¬ 
ciple requires only a CASSSCF zero-order treatment in¬ 
cluding two electrons in two MOs. We have considered 
three problems of that type. 

Bond breaking of the F2 molecule 

The F 2 molecule is a paradigmatic molecule since it 
is a case where the dynamical correlation brings a cru¬ 
cial contribution to the bonding. Despite the closed shell 
character of the wave function in the equilibrium region 
the single reference Hartree-Fock (HF) solution is un¬ 
bound (by 18 kcal/mol) with respect to the restricted 
open shell HF solution of the fluorine atoms. The 2- 
electron in 2-MO CASSCF treatment binds the molecule 
by 18 kcal/mol, but the experimental binding energy 
is much larger (39 kcal/mol). Going to a full valence 
CASSCF (14 electrons in 10 MOs) does not bring any 
improvement. The role of the dynamical correlation has 
been extensively studied and may be seen as a dynamic 
response of the lone pair electrons to the fluctuation of 
the electr ic held created by the two electrons of the a 
boncPaidO] concept of orbital breathing has been 

proposed to express the fact that the orbitals of the lone 
pairs tend to become more diffuse on the negative cen¬ 
ter and more contracted on the positive center in the 
ionic valence-bond (VB) components of the CAS. These 
dynamic relaxation processes can only take place if one 
uses non-minimal basis sets. 



FIG. 1. F 2 molecule at A=1.45 A, cc-pVDZ basis set. Con¬ 
vergence of the MRCCSD energy along the iterations. 

The calculations of F 2 were obtained in the cc-pVDZ 
basis set!^ keeping the Is electrons frozen, and accurate 
FCI estimates are taken from the work of by Bytautas et 
aP^. Figure shows an exponential convergence of the 
energy along the MRCC iterations. The here-reported 
calculation, performed in a medium size basis set, does 
not afford a sufficient hexibility to reach the experimental 
binding energy (the estimated FCI binding energy in this 
basis is Dg=28.3 kcal/mol). The potential energy curves 
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FIG. 2. Potential energy curves of the F 2 molecule, cc-pVDZ 
basis set. 



F—F distance (A) 

FIG. 3. F 2 molecule, cc-pVDZ basis set. Errors with respect 
to the FGI estimate as a function of the F—F distance. 


and the error to FCI estimate are reported, respectively 
in Figure and Figure and the estimated FCI values 
together with the error of the MRCCSD and CASSDCI 
calculations appear in Table |l] The average error is re¬ 
duced by a factor close to 6, and the NPE is only reduced 
by 40% by the MRCCSD calculations. 


The C—C bond breaking in ethane 

This calculation is performed in the 6-31G basis set, 
keeping the Is electrons frozen. The geometrical param¬ 
eters are given in Tabl e [III The potential energy curves 
are reported in Figure and the errors with respect to 
the FCI estimate appear in Figure These data show 
that the error with respect to the FCI energy is greatly 
reduced by a factor of 6 in average. According to Ta¬ 
ble [HTj the NPE goes from 2.01 mAh to 1.32 mAh for 
respectively the CASSDCI and MRCCSD approaches. 
Concerning the spectroscopic constants, the impact of 
the CC treatment is modest but goes in the right direc¬ 
tion. 


R (A) 

-E'cassdci — -E'fci 

SmRCCSD — Efci 

FCI estimate 

1.14 

19.223 

3.726 

-199.007 18 

1.20 

19.495 

3.823 

-199.048 11 

1.30 

19.825 

4.102 

-199.085 10 

1.36 

19.829 

4.045 

-199.095 17 

1.41193 

19.721 

3.920 

-199.099 20 

1.50 

19.518 

3.830 

-199.099 81 

1.60 

19.094 

3.466 

-199.095 10 

1.80 

18.038 

2.843 

-199.080 90 

2.0 

16.850 

2.034 

-199.068 82 

2.2 

16.280 

1.783 

-199.061 65 

2.40 

16.225 

1.936 

-199.058 23 

2.80 

16.055 

1.794 

-199.055 77 

8.00 

16.241 

1.893 

-199.055 45 

-^eq 

1.466 

1.465 

1.460 

^cq 

0.730 

0.739 

0.795 

Ae 

26.01 

26.91 

28.31 


TABLE I. F 2 molecule, cc-pVDZ basis set. Total energies are 
given in Ah, and the energy differences are given in mAh. 



FIG. 4. Ethane molecule, 6-31G basis set. Potential energy 
curves along the G—G bond stretching. 



FIG. 5. Ethane molecule, 6-31G basis set. Errors with respect 
to the FGI estimate as a function of the G—G distance. 
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Geometrical parameters 

C2H6 

C2H4 

C—H 

1.103 A 1.089 A 

C—C 

1.550 A 1.335 A 

H—C—C 

111.2° 

120.0° 

H—C—H 

107.6° 

120.0° 

H—C—C—H 

180.0° 

180.0° 


TABLE II. Geometries used for ethane and ethylene. 


R (A) 

£^cassdci — ^'fci 

£^mrccsd — £^fci 

FCI estimate 

4.00 

15.508 

1.834 

-79.253 166 

3.50 

15.770 

1.915 

-79.256 574 

3.00 

16.379 

2.074 

-79.268 617 

2.60 

17.037 

2.622 

-79.293 972 

2.30 

17.402 

3.009 

-79.326 999 

2.00 

17.519 

3.134 

-79.370 376 

1.90 

17.510 

3.150 

-79.385 598 

1.80 

17.482 

3.152 

-79.399 969 

1.70 

17.442 

3.106 

-79.412 107 

1.65 

17.419 

3.035 

-79.416 695 

1.60 

17.395 

3.046 

-79.419 813 

1.55 

17.371 

3.055 

-79.420 987 

1.50 

17.347 

3.062 

-79.419 613 

1.45 

17.326 

3.083 

-79.414 941 

1.40 

17.306 

3.099 

-79.406 030 

1.35 

17.291 

3.135 

-79.391 701 

1.30 

17.284 

3.153 

-79.370 480 

Ee, 

1.549 

1.550 

1.550 

^eq 

1.018 

1.017 

1.015 

Ee 

104.52 

104.99 

105.75 


TABLE III. Ethane molecule, 6-31G basis set. The FCI es¬ 
timate is the CIPSI calculation. Total energies are given in 
Eh, and the energy differences are given in mEh. 


The rotation of the ethylene molecule around its C—C 
bond 


This twisting breaks the tt bond. The calculation is 
performed in the 6-31G basis set at the geometry given 
in Table keeping the Is electrons frozen. The occu¬ 
pied MOs in the inactive space involve 10 electrons, and 
despite the modest size of the basis set one may expect 
a significant size-consistence defect of the CASSDCI re¬ 
sults, since they miss the repeatability of inactive dou¬ 
ble excitations on the SD determinants. The potential 
energy curve along the angle of rotation is reported in 
Figure and the error to the FCI estimate is reported 
in Figure From these data it appears that the global 
shape of the potential energy curve obtained using the 
CC treatment is more parallel to the FCI curve than 
using the CASSDCI approach. From Table llV) one ob- 



FIG. 6. Ethylene molecule, 6-31G basis set. Energy as a 
function of the rotation around the C—C bond. 



FIG. 7. Ethylene molecule, 6-31G basis set. Errors with re¬ 
spect to the FCI estimate as a function of the rotation around 
the C—C bond. 


serves that the error with respect to the FCI estimate is 
reduced by a factor of 6 when going from CASSDCI to 
MRCCSD. Also, the NPE is also reduced from 1.6 mEh 
to 0.3 mEh. 


Angles (degrees) 

F^cassdci — E^ci 

£^mrccsd — F^FCI 

FCI estimate 

0 

13.255 

2.935 

-78.216 340 

10 

13.132 

2.935 

-78.214 241 

20 

13.196 

2.938 

-78.208 391 

30 

13.331 

2.955 

-78.198 732 

40 

13.513 

2.991 

-78.185 373 

50 

13.750 

3.035 

-78.168 619 

60 

14.043 

3.120 

-78.149 094 

70 

14.368 

3.212 

-78.128 205 

80 

14.631 

3.258 

-78.109 498 

85 

14.694 

3.237 

-78.103 326 

90 

14.605 

3.227 

-78.100 966 


TABLE IV. Rotation of the ethylene molecule, 6-31G basis 
set. The FCI estimate is the CIPSI calculation. Total energies 
are given in Eh, and the energy differences are given in mEh. 














































































B. Two-bond breakings 


R (A) -Ecassdci — ^-FCi ^MRCCSD — Efci FCI estimate 


Three systems have been treated using a CAS with 
four electrons in 4 active MOs. Two of them concern the 
simultaneous breaking of two bonds. 

Breaking of the C=C double bond of ethylene 



FIG. 8. Potential energy cnrves of the ethylene molecule, 6- 
31G basis set. 


1.20 

9.962 

2.635 

-78.179 508 

1.25 

9.960 

2.614 

-78.200 805 

1.30 

9.961 

2.578 

-78.212 507 

1.35 

9.971 

2.578 

-78.216 869 

1.40 

9.985 

2.549 

-78.215 666 

1.45 

10.006 

2.558 

-78.210 313 

1.50 

10.028 

2.539 

-78.201 918 

1.55 

10.055 

2.517 

-78.191 365 

1.60 

10.085 

2.507 

-78.179 345 

1.65 

10.121 

2.500 

-78.166 404 

1.70 

10.158 

2.498 

-78.152 957 

1.80 

10.238 

2.528 

-78.125 775 

1.90 

10.313 

2.545 

-78.099 562 

2.00 

10.368 

2.635 

-78.075 294 

2.30 

10.310 

2.411 

-78.017 944 

2.60 

9.890 

2.228 

-77.983 952 

3.00 

9.358 

1.863 

-77.964 220 

3.50 

8.321 

1.011 

-77.956 421 

4.00 

8.616 

1.237 

-77.955 127 

-^cq 

1.362 

1.362 

1.362 

^eq 

2.043 

2.039 

2.042 

Ee 

163.48 

163.54 

164.47 



FIG. 9. Potential energy curves of the ethylene molecule, 6- 
31G basis set. Errors with respect to the FGI estimate as a 
function of the C—C distance. 


TABLE V. Dissociation of the ethylene molecule, 6-31G basis 
set. The EGI estimate is the GIPSI calculation. Total energies 
are given in Eh, and the energy differences are given in mEh. 


set at hve different geometries obtained from the equilib¬ 
rium geometry (i?e=l. 84345 A, and Zroh = 110.6°), in 
order to compare with the values of the literature.l^ The 
results appear in Table VI The benefit of the MRCCSD 
with respect to the CASSDCI treatment is significant : 
the maximum error is 1.4 mEh, better than the 6.4 mEh 
given by the Mk-MRCC treatment. This improvement 
may be due to the here-proposed treatment of the am¬ 
plitudes responsible for potential divergences. The NPE 
goes from 2 mEh to 0.7 mEh when the CC treatment is 
applied. 


The dissociation of the ethylene molecule by breaking 
the double bond was studied in the 6-31G basis set, with 
the geometry given in Table |lT| We report the potential 
energy curves in Figure and the error with respect to 
the FCI estimate in Figure The corresponding values 
appear in Table |V] Again, the error to estimated FCI 
energy is reduced by a factor of 4, but the NPF is reduced 
only by 20% with the CC treatment. 

Two-bond breaking in H2O 

This is a rather well known test problem for MRCC 
methods. The calculation is done with the cc-pVDZ basis 


R (A) 

-Ecassdci — Epci 

-E^Mk — EfCI 

E^mrccsd — E^fci 

FCI 

1 Re 

4.923 

2.909 

1.407 

-76.241860 

1.5 Re 

4.674 

4.817 

1.248 

-76.072348 

2.0 Re 

3.665 

6.485 

0.855 

-75.951665 

2.5 Re 

3.097 

5.672 

0.763 

-75.917991 

3.0 Re 

2.959 

3.987 

0.845 

-75.911946 


TABLE VI. Symmetric dissociation of the water molecule, cc- 
pVDZ basis set. The EGI total energji® is given in Eh, and 
the deviations to this reference are given in mEh. Comparison 
with Mukherjee’s state specific MRCC values (EMk — Epci) 
obtained from Ref.^^. 
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V. PROPERTIES 


VI. PROSPECTS 


A. Internal decontraction 


The method is internally decontracted. The coeffi¬ 
cients of the references as well as those of the Singles 
and Doubles change along the iterations. If the refer¬ 
ence space is a valence CAS, treating the non-dynamical 
correlation effects, the method takes care of the impact 
of the dynamical correlation on the non-dynamical part. 
The phenomenon is especially important in magnetic sys¬ 
tems where the dynamical charge polarization effects in¬ 
crease dramatically the weight of the ionic Valence Bond 
components, diminishing severely the effective energy of 
these components.!^ This effect is already present in the 
CASSDCI calculation but the MRCCSD treatment elim¬ 
inates the size consistency defect and slightly improves 
the quality of the projection of the wave function on the 
CAS. 


B. Size consistence 


The method does not introduce any unlinked diagram, 
and is therefore size-consistent. A proof of strict sep¬ 
arability has been given in the original presentation of 
the method^!. It requires that in the splitting into two 
subsystems A and B the active and inactive MOs are 
localized on one of the two subsystems A or B. Actu¬ 
ally, as occurs for the Mk-MRCC formalism, the method 
is not invariant with respect to the unitary transforma¬ 
tion of the MOs in their class (inactive occupied, active, 
inactive virtual). This dependence will be studied in a 
future work, but the error to FCI being small we do not 
expect a strong dependence on the MO definition. As 
was shown in the study of bond breakings the asymp¬ 
totic size-consistency error (which is demonstrated to be 
zero when localized MOs are used) is negligible in a basis 
of symmetry-adapted MOs. 


C. Eigenfunction of 


The here-proposed method does not provide an eigen¬ 
function of as we consider only the determinants that 
are connected by an application of H to the determinants 
belonging to the CAS. This treatment does not include 
higher excitations which will generate the full space asso¬ 
ciated with a given space part. Along all the performed 
calculations on singlet states, the order of magnitude of 
the expectation value of calculated on both the CASS¬ 
DCI and the projected MRCCSD wave functions never 
exceeded 10“^. A future work will present a solution 
working with the same restricted space but providing a 
strict eigenfunction of 


The formalism presented here allows us to conceive 
two main types of extensions for further work. The first 
one concerns the reduction of the computational cost of 
the method through various approximations, in order to 
target more realistic systems. From a methodological 
point of view, a refined treatment of the excited states 
deserves to be considered. 


A. Computational cost 

The method is extremely flexible, either on the choice 
of the reference space and/or of the excitations from it. 
One may partition the Singles and Doubles in terms of ex¬ 
citation classes. For instance the most numerous purely 
inactive excitations (2 holes, 2 particles) can be treated 
in a contracted manner, leading to a diagonal shift of 
the CASSDCI. Another possibility consists in omitting 
this class of excitation which does not contribute signif¬ 
icantly to the vertical energy differences, as exploited in 
the DDCI frameworlP^. Then, one may exponentialize 
the semi active excitations and make the DDCI method 
size consistent. As the theory is determinant based, one 
can take advantage of this flexibility to realize a CIPSI 
like selection of the dominant contributions of both the 
references and the single and double excitations. Further 
works will investigate the various possibilities such as the 
combination of MRCC with perturbation theory. 


B. Excited states 


The method is applicable to excited states using sev¬ 
eral approaches. The formalism being state specific, the 
dressing technique of the Cl matrix can be applied to any 
state dominated by the reference determinants, as long as 
a state following procedure is applied. For states belong¬ 
ing to the same symmetry, the resulting eigenvectors will 
not be strictly orthogonal but might be orthogonalized 
a posteriori. Another possibility consists in a state av¬ 
erage procedure where the amplitudes are obtained from 
the values of the quantities Xi averaged over all desired 
eigenstates: 


Em iCf 


(34) 


If one refers to the perturbative expression of the first 
order coefficients. 


A™ - 


(35) 


this approximation should be relevant when the states 
are close enough in energy. 

A recent papei!^ has proposed a generalization of this 
approach to the simultaneous treatment of several states 
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of the same symmetry. The basic ideas are the same, 
except for the fact that the extraction of the amplitudes 
is more complex. The method requires to partition the 
reference space into a main and an intermediate model 
spaces, in the spirit of the intermediate Hamiltonian for¬ 
malism. This proposal will be tested in a further work. 

VII. CONCLUSION 

This work shows the relevance of a solution previously 
proposed to the problem of the multiple parentage faced 
by all Multi-Reference treatments, as soon as the num¬ 
ber of targeted vectors is lower than the number of Ref¬ 
erences. The proposed MRCC algorithm is simple. It 
only introduces two-body excitation operators and the 
number of amplitudes to be determined is reduced to 
the very minimum. It proceeds through an iterative 
dressing of the MRCI matrix, formulated in terms of a 
standard eigenvalue equation. It is parallelizable in the 
most expensive step (the generation of the coefficients 
of the Triples and Quadruples). It is entirely decon- 
tracted and may be applied to excited states. For the 
list of benchmark studies we have performed, the results 
are extremely encouraging. The present version is state- 
specific but the principles of extension to a multi-root 
version have been formulated. This work actually opens 
into several directions, which have to be explored in the 
future. The reduction of the computational cost might be 
done using several approximations involving the selection 
of the references and/or Singles and Doubles according 
to various criteria. Furthermore, the excited states can 
be treated using different approaches, all of them being 
compatible with the here-proposed formalism. 

On a different perspective, the multiple parentage 
problem, which was faced here in the purpose of building 
a logically consistent computational tool to go in the di¬ 
rection of the exact solution, also concerns the building of 
rational valence-only effective Hamiltonians. In such an 
approach, the idea is to map the information coming from 
a sophisticated treatment into a minimal effective Hamil¬ 
tonian, the parameters of which should be as physically 
meaningful as possible.l^ We believe that the solution 
we proposed to the multiple parentage problem offers a 
rational solution to this reduction of information. This 
remark illustrates the intrinsic link between the two main 
tasks of Quantum Chemistry, namely the production of 
physically grounded interpretative models on one hand 
and the conception of rigorous computational tools. 
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